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Nonlinear dynamics of the ion Weibel-filamentation instability: an analytical model 
for the evolution of the plasma and spectral properties 
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We present a predictive model of the nonlinear phase of the Weibel instability induced by two 
symmetric, counter-streaming ion beams in the non-relativistic regime. This self-consistent model 
combines the quasilinear kinetic theory of Davidson et al [Phys. Fluids 15, 317 (1972)] with a 
simple description of current filament coalescence. It allows us to follow the evolution of the ion 
parameters up to a stage close to complete isotropization, and is thus of prime interest to understand 
the dynamics of collisionless shock formation. Its predictions are supported by 2-D and 3-D particle- 
in-cell simulations of the ion Weibel instability. The derived approximate analytical solutions reveal 
the various dependencies of the ion relaxation to isotropy. In particular, it is found that the influence 
of the electron screening can affect the results of simulations using an unphysical electron mass. 

PACS numbers: 


First-principles kinetic simulations of plasma collisions 
governed by electromagnetic effects are now made possi¬ 
ble using massively parallel particle-in-cell (PIC) codes, 
hence paving the way to quantitative modeling of a num¬ 
ber of high-energy astrophysical scenarios p], Q . The tur¬ 
bulent shocks possibly arising from plasma instabilities in 
these systems are believed to be responsible for the gener¬ 
ation of nonthermal particles and radiation 00. In this 
context, many simulation studies have demonstrated the 
ability of the Weibel-filamentation instability to 

provide the electromagnetic turbulence required for effi¬ 
cient dissipation of the flow energy and Fermi-type accel¬ 
eration processes [2, [l2|, 13]. These numerical advances 
go along with experimental progress towards the laser- 
driven generation of collisionless turbulent shocks in the 
laboratory [3fE3. 

Collisionless shocks developing in electron-ion plasmas 
may be of laminar or turbulent nature depending on the 
type (electrostatic or electromagnetic) of the dominant 
underlying instability jl8]. In this work, we concentrate 
on initially unmagnetized electron-ion systems whose col¬ 
lective dynamics is eventually ruled by the electromag¬ 
netic ion Weibel instability, which may evolve into a tur¬ 
bulent shock. While thisproblem has inspired a number 
of numerical studies [l|, [2, Qji, fl6l4l9l ] , there is as yet no 
analytical model of the nonlinear evolution of the ion 
Weibel instability leading to shock formation. Our goal 
is to provide such a description within the simplifying as¬ 
sumption of homogeneous and infinite colliding plasmas 
of equal densities and temperatures. Our paper is orga¬ 
nized as follows. In Sec. HI we first analyze the results 
of a reference PIC simulation, pointing out the transi¬ 
tion from the early-time electron-driven phase, associ¬ 
ated with various fast-growing modes, to the ion-driven 
phase ruled by the ion Weibel instability. In Sec. [Ill 
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we present a set of quasilinear equations describing the 
evolution of the ion parameters in the time-varying mag¬ 
netic turbulence generated by the Weibel instability J20j . 
Approximate relations between the plasma and spectral 
parameters are obtained and successfully confronted to 
a number of PIC simulations. Our model is then made 
self-consistent by the inclusion of a simple description of 
current filament coalescence. The derived analytical so¬ 
lutions are shown in good agreement with PIC simulation 
results. Our concluding remarks are given in Sec. IIII1 

I. TRANSITION BETWEEN ELECTRON AND 
ION INSTABILITY REGIMES 

The PIC simulations presented in this study have been 
performed using the parallel code calder [21|. All of 
them consider plasma systems made of two symmetric, 
counter-streaming electron-ion beams of uniform den¬ 
sity and velocity profiles. Both the electrons and ions 
initially obey Maxwell-Jiittner distributions of temper¬ 
atures T e = Ti and drift velocities v e = Vi = zb^ox. 
Periodic boundary conditions are used for the fields and 
particles in all directions. Our reference two-dimensional 
(2-D) simulation, examined below, is characterized by 
T e ,i/Tn e c 2 = 0.01 and Vo/c = 0.2. We use singly-charged 
ions with a reduced mass rrii/m e = 100 to alleviate 
the computational cost. The density of each beam is 
normalized to the total electron density (n e ^ = 1/2). 
The mesh sizes are Ax = Ay = 0.1 c/uo pe , where uo pe is 
the total electron plasma frequency. The time step is 
At = 0.069^“^. The domain size is 1024Ax x 1024 Ay, 
that is, 102 Ac/uj pe x 102 Ac/uj pe . Each cell initially con¬ 
tains 50 macro-particles per species, yielding a total num¬ 
ber of about 2 x 10 8 macro-particles. Third-order weight 
factors are employed along with Esirkepov’s current de¬ 
position scheme |22l ]. 

In the relativistic, cold electron-positron systems con¬ 
sidered in Refs. [23|, [24[, the Weibel filament at ion was 
found to prevail from early on. By contrast, the non- 
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FIG. 1: (a) Growth rate Y /uj pe in the (k x , fej-plane for the plasma parameters rrn = 100?7ie, v e #§ Vi = ±0.2c and T e — Ti — 
m e (?/ 100. The real frequency vanishes for all the illustrated unstable modes, (b) Temporal evolution of the space-integrated 
E x (black), E y (red) and B z (blue) electromagnetic energies, (c) Total ion density (normalized to the total initial density n- 0 ^) 
at ujpet = 690. (d) p x — p y electron phase space at c o pe t — 690. (e) Magnetic spectrum, e\B z (k x , k y )\/m e ^pe (in log 10 scale), at 
ujpet = 690. (f) Transverse lineouts of the E y and B z helds and of the particles’ mean x-velocity at uj pe t — 690. 


or weakly-relativistic, warm electron-ion systems treated 
henceforth usually experience a progressive transition 
from an early phase ruled by electron-driven modes to a 
regime ruled by the slower ion Weibel instability. Our ref¬ 
erence simulation illustrates this transition. Figure QJa) 
displays the theoretical growth rate map in the (k x ,k y ) 
plane, computed from the electromagnetic dispersion re¬ 
lation (see the Appendix) with the initial plasma pa¬ 
rameters. The largest growth rate, T max /u) pe ~ 0.18, 
corresponds to a longitudinal electrostatic mode located 
at ( k x ,k y ) ~ (2.5, 0)uj pe /c. The electromagnetic Weibel 
modes are located around the k y - axis, with a maximum 
growth rate Y/uj pe ~ 0.1 reached at the purely trans¬ 
verse wave vector (fc x , k y ) ~ (0, l)w pe /c [Fig. Oja)]. Note 
that the growth rate map’s fastest-growing mode associ¬ 
ated with a given wave vector. In the present case, these 
dominant modes are essentially driven by the electrons 
and, owing to the symmetry of the system, present a van¬ 
ishing real frequency. We verified that the calculation 
performed with immobile ions give negligible differences 
in the growth rates. 

The predicted electrostatic character of the early in¬ 
teraction phase is confirmed by the time history of the 
integrated electromagnetic energies [Fig. [H[b)] . At early 
times ( uj pe t < 40), the E x energy is the dominant one, 
growing exponentially at an effective rate Y/uj pe ~ 0.14 
consistent with linear theory. The associated longitudi¬ 
nal instability saturates at uj pe t ~ 40 and subsequently 
decays away [Fig. |T](b)] . The magnetic B z energy then 
takes over, growing at a rate Y/uj pe ~ 0.08 close to the 
fastest Weibel mode of Fig. Oja). A first magnetic sat¬ 


uration occurs at uo pe t ~ 80, which marks the end of 
the electron-governed Weibel instability. The instability 
then switches to an ion-driven regime [iima ini¬ 
tially characterized by exponentially increasing B z and 
E y energies (while the E x energy keeps stagnating at 
a low level). This growing phase comes to an end at 
u: pe t ~ 500, at which time the instability enters its non¬ 
linear saturation phase, further studied in the following. 
This phase exhibits stagnating field energies but also, as 
analyzed below, continuously evolving plasma and spec¬ 
tral distributions. 

The transverse character of the magnetic instability 
is evidenced in Fig. DJc), which displays (in log 10 scale) 
the magnetic spectrum \B z (k x ,k y )\ at uj pe t = 690. The 
spectral energy is concentrated along the k y - axis, peak¬ 
ing around (. k x ,k y ) ~ (0, 0.7)u pe /c. The dominant ion 
contribution to the instability is demonstrated by the 
transverse field and current lineouts plotted in Fig. [IJd). 
The total ion current appears more strongly modulated 
than the electron current, so that the magnetic fluctua¬ 
tions are mainly induced by the ions. This feature, shown 
here at uj pe t = 690, is found to hold from uj pe t ~ 300 
onwards. As expected in its nonlinear regime, the ion 
Weibel instability generates, besides current modula¬ 
tions, significant ion density fluctuations in the trans¬ 
verse direction [Fig. [U(e)]. These ion density filaments 
also present longer-wavelength longitudinal modulations 
that can be ascribed to magnetic coalescence processes 
mill. The latter will be shown to govern the non¬ 
linear evolution of the magnetic turbulence. The partial 
electron screening of the ion density filaments accounts 
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for the correlated growth of the E y and B z energies seen 
in Fig. DJ(b). From the lineouts of Fig. QJd), the E y fluc¬ 
tuations have a typical wavelength twice smaller than the 
B z fluctuations, which suggests that an approximate bal¬ 
ance between the transverse electric and magnetic forces 
is established on the electrons [27j. From this reason¬ 
ing, the transverse electric field is expected to scale as 
E y ~ 2 E~d y Al, where A x is the vector potential. As¬ 
suming A x ~ Bo cos(k y y)/k y (where Bo is the magnetic 
field amplitude and k y the dominant transverse wave 
number), one predicts the transverse electric field am¬ 
plitude Eo ~ 2m e k F>q. Using eBo/m e uj pe = 0.15 and 
k y c/co pe = 0.7, one finds eEo/m e cuj pe = 0.016, in fair 
agreement with Fig. Hid). 

Due to the ensemble of electron-driven instabilities de¬ 
veloping at early times [Fig. QJa)], the electrons turn out 
to be essentially isotropized in the ion Weibel regime, 
as shown by their p x — p y phase space at uj pe t = 690 
[Fig. [Hf)]. More quantitatively, the average x- velocity 
of each electron beam has dropped from \v e \ = 0.2c to 
\v e \ ~ 0.06c at the beginning of the nonlinear ion Weibel 
phase (uj pe t ~ 500), while its temperature has increased 
to T e /m e c 2 ~ 0.04 [see Figs. [5][d,f)]. 


II. QUASILINEAR MODEL OF THE ION 
WEIBEL FILAMENTATION 

A. Dominant unstable mode 

A major assumption of our model is that the main 
properties of the magnetic spectrum can be related to the 
instantaneous plasma parameters. In order to derive such 
a relation, let us first examine the linear characteristics of 
the purely transverse ion Weibel instability. Its general 
dispersion relation follows from taking 0 = 7r/2 in Eq. 

(ED): 


(^-klc^Syy-U 4 4 = 0 . ( 1 ) 

In the case of symmetric counter-streaming plasmas, the 
off-diagonal tensor element vanishes .yielding the well- 
known simplified dispersion relation p] 


imaginary character of the Weibel modes (cu = iF with 
r > 0), the dispersion relation can be recast as [20] 


k y c 2 + T 2 + ^cj 2 s - ^u; 2 s (a s + 1)5R [1 + i s Z (£ s )] — 0 . 

s s 

(4) 

In the above equation, Z denotes the plasma disper¬ 
sion function (28j, uo ps the plasma frequency of the sth 
species, £ s = iy/m s /2T ys F/k y , 5R the real part and a s the 
anisotropy ratio of the sth species, defined as 


a s 


TYlsV 'o + T s: 


- 1. 


L sy 


(5) 


We will now assume that, independently of their ini¬ 
tial distribution, the electrons are almost completely 
isotropized (a e ~ 1) in the nonlinear ion Weibel phase. 
Furthermore, they will be assumed hot enough so that 
|£ e | 1, allowing us to use the small-argument expan¬ 

sion of Z\ 


z{£) = *v / 7T|^|exp(-^ 2 ) - 2£ + 0(£ 3 ). (6) 


To leading order, we obtain \£ e Z(£ e )\ ~ \\Zr^ e \ <C 1. 
For instance, solving Eq. with the plasma parame¬ 
ters measured in the reference simulation at uj pe t = 500 
( T ey 22 Ti y 22 0.04m e c 2 leads to |£ e | ~ 0.1 and \£iZ\ ~ 
0.15. This term will thus be neglected in the bracketed 
factor of Eq. (J4]) . 

The ion response will be assumed to fulfill |&| < 1, 
so that we retain the leading term of &>£(&) ~ y/i in 
Eq. (J2J). There results the approximate dispersion rela¬ 
tion 


+ 1) 


I fe 2 

ox 1 . 1 + s 

Z'-L iy 


u2 W pe«e 

S r 2 r 2 


The growth rate is readily solved as 

-.2 22 i2 22 


r ~ 


9 rp L.2 2 _ jL2 '2 

^J-iy | L | ^max^ 

\Ky\- 


Trrrii 


uli( a i +1) 


= 0. 

(7) 

(8) 


where & max denotes the upper bound of the Weibel- 
unstable domain: 


u 2 £ xx - kyC 2 = 0 . (2) 

For the sake of simplicity, we will restrict our analy¬ 
sis to non-relativistic systems initially described by two- 
temperature, drifting Maxwellians of the form 


m 8 (v x -v 3 ) 2 m s v y 

2 tTx ~ZEy J ’ 

(3) 

where m s is the mass, v s is the drift velocity and T xs 
and T sx are, respectively, the longitudinal and transverse 
temperatures of the 5th plasma species. Using the ex¬ 
pressions (1A3D . (1A4D and (IA5D . and exploiting the purely 


/i 0 ) (v) = 


m s 


27 t^/T sx T s . 


: exp 


^max — C -\JUJp e CL e H - UJ p ^CLi . (9) 

Equation (|5]) is formally similar to that derived for the 
electron Weibel instability in the weak growth rate limit 
fjoj . The only difference lies here in /c max , which involves 
both the ion and the electron anisotropies. 

The accuracy of the estimate, Eq. (j8]), is illustrated 
in Fig. [2] for the plasma parameters Vi = ±0.2c, v e = 
0, T ex ^ y = Ti X:y = 0.05m e c 2 and mi = 25m e . Correct 
agreement is found between the exact and approximate 
growth rate curves, with respect to both their general 
shapes and the location and amplitude of the dominant 
mode. 
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FIG. 2: ^-dependence of the Weibel growth rate for the 
plasma parameters Vi = =L0.2c, v e = 0, T ex , y = Ti X , y = 
0.05 m e c 2 and rrii/m e = 25: comparison between the exact 
solution of Eq. 0 (solid line) and the estimate, Eq. Q. 


In the following, we will assume that most of the mag¬ 
netic energy is contained in the spectral region surround¬ 
ing the fastest-growing wave number k sat . Making use of 
Eq. (j8j), the solution of d ky Y = 0 is 

7 ^max / in \ 

ksat ' ( 10 ) 

To further simplify our analysis, we will henceforth use 
the approximation fc sat — 0.5M max . Moreover, we will 
take the large-ion-anisotropy limit oq 1, consistently 
with our focus on the isotropization process of initially 
highly anisotropic ion populations. The oq 1 limit 
will be assumed to hold (at least marginally) up to a 
stage close to ion isotropization (a* >2). In an inhomo¬ 
geneous system, the latter stage should approach shock 
formation. There follow the estimates 


kinetic theory 

dt{fs)(v,t) = ~iJ2 


pi I 


3 k 

K'V 


lion s m s c 2 kl 


k 

X [ kyV X 0 Vy T (i^-ky T kyVy^dv 

kyV X 3y y T (iT ky kyVy'jdy 


X 


^ ky ky v y 


(fs)(v,t), (14) 


where fi o is the magnetic permittivity of vacuum. We 
recall that the quasilinear kinetic theory for the Weibel 
instability is valid provided \£ s \ < 1, as already assumed 
in the previous section. The ability of the quasilinear the¬ 
ory to capture the nonlinear evolution of the Weibel in¬ 
stability has been demonstrated by Davidson et ol. [20j . 

Assuming that the ion distribution functions remain of 
the bi-Maxwellian form, Eq. (J3j) , the three first moments 
of Eq. (fTD) give a set of differential equations on the 
mean ion drift velocities and temperatures: 


S^t^S — ^ ^ 


Wps Vs ft [1+ £,£(&)] 9t|gfc « |2 , (15) 


_ k 2 r 2 T 

fc„ K v C lsy 


s d t T sy = V + 1)3? [1 + £ a Z(i .. 

K y c n o 


Mo 
dt\B k J 2 


(16) 


s d t K sx = -J2 ^|3? [2(a s + 1)(1 + 6^(6)) - 1] 

ky y c 

dt\B ky \ 2 


Mo 


(17) 


k 


max 


WpiV^/^ i 


r 


ky 





(11) 

( 12 ) 


with K sx = T sx + m s v 2 , the x -momentum flux. We have 
also exploited the relation 

d t \B ky \ 2 = 2Y ky \B ky \ 2 . (18) 


(13) 


Approximate solutions of the quasilinear equations 


Therefore, for ai 1, it is found that & depends only 
on the ratio (k y /k mdiX _) 2 . Anticipating on the next sec¬ 
tion, we define £ sat = £i(k y = fc sat ) ~ 0.4i. The 
small-argument expansion of Z{fii) in Eq. (J4]) is then 
marginally valid. 


B. Temporal evolution of the plasma parameters 


In order to make analytical progress, we will make use 
of the approximation 


Edr^ 1 +&.*&,)] 


dt\B ky \ 2 


mi 


m+UtZ(Ut)]S p , 


(19) 


1. Quasilinear equations 

The above formulae will serve to relate to the dominant 
wave vector, fc sat , to the ion anisotropy ratio deduced 
from the spatially averaged ion distribution functions, 
(/s)(v, t). The evolution of the latter due to non-resonant 
wave-particle interaction in the Weibel magnetic turbu¬ 
lence will be described in the framework of quasilinear 


where Z{ is the ion charge number and the spectral pa¬ 
rameter S p is defined as 

5 p= e 2 E%^ = e2 E I 4 U 2 , ( 2 o) 

ky y ky 

which is homogeneous to the square of a momentum 
and where the sum runs over the positive and negatve 
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FIG. 3: Temporal evolution of the quantity Oi defined by 
Eq. (l24|) in the reference simulation (v% — v e — ±0.2c and 
nrii — 100m e ). 


For practical reasons, we define the parameter 

Ke^OiT^ +T^ +m iV f )2 , (28) 

which has the dimension of an energy, and where = 
vq. For the typical value Oi = 2, we have K 2 = 
2± ± rriiV ^ 2 . The integration of Eq. ([27]) is 

straightforward, yielding 


rp _ 

± ^y ± iy 


. K 2 - 

+ T ln 


/ K2 - fljTfr \ 
\K 2 



rrii 


C s p -s f). 


(29) 

The Taylor expansion of the logarithmic term for 
T iy /K 2 < 1/di 1 (high-anisotropy limit) gives, to lead¬ 
ing order, 


wavevectors. Within the range of validity of Eq. m, 
the factor Oii = 1 ± £ S at^(£sat) is a constant (oti ~ 0.5), 
independent of the plasma parameters. Equations (fT5]h 
(usd and CZD can then be recast as 


d t Vi = -—ai^-d t S p , 
mi T iy 

Z 2 

d t T iy = — —QLi^ai ± 1 )dfS p , 

TTli 

z 2 

dJCx =-i- (2a iffli + 2a, - 1) d t S p . 


( 21 ) 

( 22 ) 

(23) 


An additional simplification can be made upon real¬ 
izing that the ratio of Eqs. ([22]) and ([23]) is essentially 
constant in the a* 1 limit: 


dfKi x _ q ± 1 

dtTiy OLi ((Li ± 1 ) 


- -2 



(24) 

(25) 


The validity of this approximation is supported by Fig. [3j 
which plots the time history of Oi from the reference sim¬ 
ulation. This quantity is seen to rapidly drop down to a 
saturated value ~ 2 once the ion Weibel instability en¬ 
ters its nonlinear phase ( uj pe t > 500). At later times, Oi 
slowly decreases (reaching ~ 1.8 at uj pe t = 2700) as a 
result of decreasing ion anisotropy [see Fig. E3[e)]. We 
will then assume that 0i is a quasi-constant (to leading 
order in 1/a*) and neglect its time derivatives. Equation 
(|24|) can then be readily integrated, giving 


Ti V ~ \ T.;/ ± 2 


( 0)2 . 


iy 


’—aiK 2 (S p - 4 0) ) . (30) 

rrij \ J 


Inserting Eq. ®HD into m yields 


, Zf 

OtVi ~- -oti 

rrii 


^i Of Sp 


T% )2 + 2: 


(Sp - 4° } ) 

which may be readily integrated as 


(31) 


Vi ~ v 0 exp 


-2 


T iy (S p )-T^ 

K 2 


(32) 


where T iy verifies Eq. (TITfl) . Combining Eqs. (T2TTT) and 
m allows us to solve for Ki X : 


K lx ^K 2 -ejT^ 2 + 2?ta i K 2 ( 


There follows the anisotropy ratio 

K 2 


q _ o(°) 

o p Op 


- 2 . 


). (33) 


T^ 2 + 2^-a i K 2 {S p -S^ 


(s P - 4 0) ) 


(34) 


and the spectral parameter 


Q _ C(0) ^ 

P ^ “2 ZfoiK 2 


K 2 

2 + (ij 


- T, 


( 0)2 


(35) 


K ix = K^-e i (r iy -T^) , ( 26 ) 


where the notation stands for X(t = 0). Plugging 
Eq. ([26]) into ([22]) leads to 


L iy 


/r(°) — 0 (T- -T ( 

^ix U i\- L iy ± i 


(0) 

iy 


-d t Ti 


t-Liy 


Z ? 


rrii 


ctidtSp 


(27) 


To summarize, using non-resonant quasilinear theory 
[20| . we have expressed the ion parameters Vi, cii Ki x 
and in terms of the instantaneous spectral parameter 
independently of the time history of the ion Weibel- 
governed system. The reader should be reminded that 
these approximate relations are valid in the limit of es¬ 
sentially isotropic electrons and highly anisotropic ions. 
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(a) Electromagnetic energies 


(b) Spectral parameter S p 




(d) Mean ion ^-velocity (e) Ion anisotropy ratio 




(c) Dominant wavevector 



(f) Transverse temperatures 



FIG. 4: 2-D simulation with initial conditions: v% = ±0.2c, v e = 0, T e ^/m e c 2 = 0.01, mi/m e a* 100. (a) Temporal evolution 
of the mean electromagnetic energies (normalized to m e c 2 u pe /c 2 ). (b) Temporal evolution of S p /(m e c) 2 from the simulation 
(black solid line), (c) Temporal evolution of k sa tc/uj pe maximizing the l±-spectrum from the simulation (black solid line) and 
from Eq. (1101) (black dashed line). Temporal evolutions of (d) Vi/c, (e) ai and (f) Ti, ey /m e c 2 from the simulation (black solid 
lines) and from Eqs. (l32l) . (l34l) and (l30l) . respectively (black dashed lines), extracting S p from the simulation [black solid line 
of (a)]. The analytical predictions, Eqs. l|521) . iP |l - l|5Tj> . are superposed as red dashed lines. 


C. Comparison with PIC simulation results 

This section gathers the PIC simulation results ob¬ 
tained using a variety of periodic (2-D or 3-D) geometries 
and plasma parameters, and confronts them to the above 
analytical expressions. The latter will be computed using 
the PIC-predicted values of either the spectral parame¬ 
ter, S p (t), or the ion anisotropy ratio, a*(t). 


1. 2-D periodic simulations 


For the 2-D PIC simulations, the spectral parameter 
S p is computed using a discrete Fourier transform in the 
transverse y direction, averaged along the x direction: 


n2D 

Op 


2^(\DFT y (B z )\ 2 ) x 

Z_—/ b>2 


(36) 


Since we neglect the electron anisotropy in our quasilin- 
ear model, let us first consider a plasma system with 
initially isotropic electrons. This configuration is ex¬ 
emplified in Figs. Ha-f), which gather the results of a 
simulation run with Vi/c = ±0.2, v e = 0, T e ^/m e c 2 = 
0.01 and rrii/m e = 100 (black solid curves). The nu¬ 
merical resolution of the corresponding dispersion rela¬ 
tion (not shown) predicts that the fastest-growing mode 
(rmax /^pe — 0.013) is of the Weibel kind, and that there 
is no unstable longitudinal mode. These predictions are 


confirmed by Fig. |U[a), which shows that the system’s 
evolution is ruled from the start by the magnetic field 
growth. The B z energy exponentially increases during 
the time interval 200 < uj pe t < 600, at an effective rate 
r/cjp e ~ 0.01, close to the theoretical value. At the sat¬ 
uration time ( uj pe t ~ 600), the B z energy is about 40 
times larger than the E y energy. Later on, although the 
field energies vary very weakly, the ion parameters (vi, 
Gq, Tiy) continuously evolve [Figs. IU[d-f)]. These distinct 
behaviors point to a time-changing magnetic spectrum, 
as confirmed by the increasing (resp. decreasing) trends 
of S p (resp. fc sat ) observed in Figs. [H(b,c). 

Let us confront these results to our quasilinear model. 
Its validity is first verified by solving the dispersion rela¬ 
tion with the plasma parameters measured at the satura¬ 
tion time. We obtain |£ sat | ~ 0.5 < 1, in good agreement 
with the approximate value of 0.4 obtained in Sec. Ill Al 
This should be contrasted with the value £ sat ~ 1.2 > 1 
associated with the initial state of the system. The 
dashed black lines in Figs. |U[d,e,f) plot Vi(t), cq(£) and 
T iy (t) as predicted by Eqs. (I32]h (tMl) and (tHUl) using S p (t) 
from the simulation [Fig. ffl[b)]. Satisfactory agreement 
with the PIC results is found for the three curves. The 
weak variations of vi are well reproduced over the whole 
simulation time. More interestingly, the pronounced vari¬ 
ations (by about an order of magnitude) of and Ti y are 
quantitatively captured up to tuj pe ~ 1400. Later on, our 
model underestimates by a factor of ^ 2 the increase in 
Ti y , and consequently, overestimates a* by the same fac- 
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(a) Electromagnetic energies (b) Spectral parameter S p 






(c) Dominant wavevector 



(f) Transverse temperatures 



tCO ne 


FIG. 5: 2-D simulation with initial conditions: Vi — d=0.2c, v e = ±0.2c, T e ^/m e c 2 = 0.01, rrii/m e = 100. (a) Temporal 
evolution of the mean electromagnetic energies (normalized to m e c 2 u: 2 e /c 2 ). (b) Temporal evolution of S p /(m e c ) 2 from the 
simulation (black solid line), (c) Temporal evolution of k sa tc/uj p e maximizing the B z -spectrum from the simulation (black solid 
line) and from Eq. (HOl) (black dashed line). Temporal evolutions of (d) Vi/c, (e) a* and (f) T e ^y/m e c 2 from the simulation 
(black solid lines) and from Eqs. (l32l) . (l34|) and (l30l) . respectively (black dashed lines), extracting S p from the simulation [black 
solid line of (a)]. The analytical predictions, Eqs. (Pl). pi)-p|). are superposed as red dashed lines. 


(a) Electromagnetic energies 



(d) Mean ion ^-velocity 



(b) Spectral parameter S p 



(e) Ion anisotropy ratio 



(c) Dominant wavevector 



(f) Transverse temperatures 



FIG. 6: 2-D simulation with initial conditions: Vi — ±0.4c, v e — 0, T e ^/m e c 2 = 0.01, rrii/m e = 1836 (a) Temporal evolution 
of the mean electromagnetic energies (normalized to m e c 2 uj pe /c 2 ). (b) Temporal evolution of S p /(m e c) 2 from the simulation 
(black solid line), (c) Temporal evolution of k sa tc/u pe maximizing the i^-spectrum from the simulation (black solid line) and 
from Eq. (HOl) (black dashed line). Temporal evolutions of (d) Vi/c, (e) a% and (f) T e ^y/m e c 2 from the simulation (black solid 
lines) and from Eqs. (l32l) . (l34l) and (l30l) . respectively (black dashed lines), extracting S p from the simulation [black solid line 
of (a)]. The analytical predictions, Eqs. Pi). (EMm, are superposed as red dashed lines. 


tor. The accuracy of Eq. (USD, giving fc sat as a function of nq (here taken from the PIC curve in Fig. S^e)), is 
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illustrated by the dashed black line in Fig. [4(c). Good 
agreement is observed between the PIC and approximate 
curves of fc sat (t). 


Let us now return to the reference simulation of v e = 
Vi = ±0.2c, the other parameters being identical to those 
of the previous case. The corresponding results are dis¬ 
played in Figs. [5(a-f). The strong deceleration of each 
electron beam is illustrated in [5(c). At the saturation 
time of the ion Weibel instability ( tu pe ~ 500), the elec¬ 
tron drift velocity, v e , has decreased by more than a fac¬ 
tor of 2 and the electron anisotropy ratio, a e , has dropped 
to a e ~ —0.2 (not shown). Using the instantaneous 
plasma parameters (2% ex ~ 0.01 ra e c 2 , T^ ey ~ 0.03ra e c 2 , 
Vi ~ ±0.2c and v e ~ ±0.08c), the dispersion relation 
gives |£ e | ~ 0.04 <C 1 and |^| ~ 0.4, the latter value 
closely matching the theoretical expectation. Plugging 
the simulated S p (t) values into Eqs. ([32]) . (]34l) and ([30]) 
yields approximate curves (dashed black lines) that re¬ 
produce the PIC curves (black solid lines) to within a 
factor of - 2 (for a* and Ti y ). Furthermore, the fc sat esti¬ 
mate, Eq. m, is seen to underestimate the simulation 
values by 40% for tuj pe > 500 [Fig. [5(b)]. 


Figures [6( a-f) show the results of a simulation run with 
a realistic proton mass, rrii = 1836m e , and Vi = ±0.4c 
and v e = 0. The other parameters are kept unchanged. 
A cruder discretization was employed for this simulation 
(.Ax = Ay = 0.2 c/u pe ). The theoretical growth rate 
map shown in Fig. 0 predicts that the system is initially 
dominated by a longitudinal electrostatic mode of growth 
rate r max /u;p e ~ 0.04, wave number k x c/uj pe ~ 2.8 and 
phase velocity uo/k ~ =L0.4c. This is indicative of a Bune- 
man instability, driven by a relative electron-ion veloc¬ 
ity exceeding the electron thermal velocity (v te = 0.1c). 
These predictions account for the initially dominant E x 
energy observed in Fig. Ha). Following an exponentially 
growing phase (at a rate Y/uj pe ~ 0.03, comparable to 
the theoretical value), the E x energy saturates around 
tuj pe ~ 200, and rapidly decays away. The B z mag¬ 
netic energy associated with the Weibel instability over¬ 
comes the E x energy at tuj pe ~ 400, before saturating 
at tuj pe ~ 1000. As in the previous cases, the magnetic 
energy remains approximately constant after saturation. 
The evolution of fc sat is well reproduced by Eq. m 
for tuj pe > 1000 [Fig. [6(c)]. Furthermore, the temporal 
evolutions of Vi, di and 7% in the simulation are well re¬ 
produced by the quasilinear theory over the whole simu¬ 
lation time [Figs. [6(d,e,f)]. Surprisingly good agreement 
is found during the early Weibel phase, although the 
\£i\ < 1 condition is not fulfilled at the saturation time: 
|£i| ~ 1.2 is then obtained for the measured parameters 
Tex = T iXt y = 0.2 m e c 2 , T e:iy = 0.3 m e c 2 , Vi = ±0.4c, 
v e = 0. Later on, |^| steadily decreases below unity, 
reaching |^| ~ 0.8 at tuj pe = 4000. 


FIG. 7: Growth rate T(k x , k y )/u P e of a colliding-beam system 
with Vi = d=0.4c, v e m 0, T ep /m e c 2 - 0.01 and mni/m e = 

1836. The real frequency vanishes for the unstable modes 
around k x — 0 (left panel), while the fastest-growing modes 
around (k x ,k y ) ~ (0.3, 0)c o pe /c (right panel) propagate at a 
phase velocity of — 0.4c. 


2. 3-D periodic simulations 

The above quasilinear equations can be readily gener¬ 
alized to 3-D systems, where the Weibel instability de¬ 
velops in the y — z transverse plane, given the following 
definition of S p : 

o3D _ r 2 V (\ DFT vABy)\ 2 + | DFTy, z (B z )\ 2 ) x 

D P e 2—/ 1.2 I 1.2 

ky,k z Ky + K z 

(37) 

Moreover, the dispersion relation, Eq. 0 , remains un¬ 
changed when shifting from 2-D to 3-D non-relativistic 
bi-Maxwellians. 

Figures [8(a-f) present the results of a 3-D periodic 
simulation using the same plasma parameters as in 
Figs. [8(a-f): v e = Vi = =L0.2c, T e ^ = 0.01m e c 2 and 
rrii = 100m e . The simulation domain has dimensions 
102.4c /uj pe x 96c/cj pe x 9 6c/uj pe with the discretization 
Ax = Ay = Az = 0.2 c/uj pe . Each cell is initialized with 
30 macro-particles per species. 

The 2-D and 3-D simulations give very similar results 
with respect to both the electromagnetic and kinetic 
quantities (Figs. [5] and [8]). A somewhat surprising find¬ 
ing is that the 2-D simulation predicts a slightly faster 
increase in the transverse ion temperature,! 1 ^. The ob¬ 
served overall agreement between 2-D and 3-D simula¬ 
tions is consistent with Ref. [33|, where it was demon¬ 
strated that the multidimensional physics of unstable 
two-stream systems is well captured by 2-D simulations 
resolving the drift (longitudinal) direction. Finally, as 
in the 2-D case, the predictions from quasilinear theory 
reasonably match the 3-D simulation results during the 
ion Weibel-saturation stage (tu pe > 400). 
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(a) Electromagnetic energies 



(b) Spectral parameter S p 



(c) Dominant wavevector 



• 3D PIC 

...Quasilinear theory, a.»1 

--2D PIC 

—Coalescence model 

V ?\"—... 

;~- :x 



°0 200 400 600 800 1000 1200 


(f) Transverse temperatures 





FIG. 8: 3-D simulation with initial conditions: Vi = ±0.2c, v e — 0, T e ,i/m e c 2 = 0.01, rrii/m e — 100. (a) Temporal evolution of 
the mean electromagnetic energies (normalized to m e c 2 u pe /c 3 ). We have checked that the B y and B z histories exactly coincide. 
The B z history from the 2-D simulation is superposed as a blue dotted-dashed line, (b) Temporal evolution of S p /(m e c ) 2 , 
defined by Eq. (l20l) from the simulation (circles), (c) Temporal evolution of k sa ,tc/uj pe maximizing the i^-spectrum from the 
simulation (circles) and from Eq. CGD (red dashed line). Temporal evolutions of (d) Vi/c, (e) cq and (f) T e ^ y /m e c 2 from the 
simulation (solid lines) and of Eqs. (l32l) . (l34l) and (l30l) . respectively (black dashed lines) measuring S p from the simulation 
(circles of (a)). The 2-D PIC simulation results of Fig. [5]are superposed as black dotted-dashed lines. The analytical predictions 
of Eqs. (l52l) . (I54l) - (l57l) are superposed as red dashed lines. 


D. Coalescence-driven spectral dynamics 

From the quasilinear theory of the transverse Weibel 
instability, we have derived simple analytical relations 
between the ion parameters and the spectral quantity 
S p = e 2 E k \A k \ 2 . These equations have been shown 
to match the PIC simulation results for various plasma 
parameters, provided the ion anisotropy remains large 
enough in the nonlinear stage. A closure relation rela¬ 
tion, however, must be provided to get a fully predictive 
model. 

Previous studies of the Weibel instability revealed that 
the nonlinear filament dynamics is subject to secondary 
processes, such as kink instabilities [34[, which generate 
k x 7 ^ 0 modes in the magnetic spectrum, or filament co¬ 
alescence [26|, l35l IdT^ . The latter mechanism originates 
from the partial neutralization by the background elec¬ 
trons of the ion current filaments formed in the nonlinear 
stage. This results in a nonvanishing magnetic attrac¬ 
tion between neighboring filaments of same current sign, 
which then tends to coalesce. Each merger generates a 
larger filament, of roughly twice the size of the primary 
filaments, thus leading to increasingly \ow-k y modes in 
the magnetic spectrum. The close agreement between 
the above simulations and the quasilinear theory of the 
purely transverse ion Weibel instability indicates that the 
long-term dynamics of the latter is mostly governed by 


the evolution of the k x = 0 modes, and therefore by coa¬ 
lescence effects. In the framework of our model, this pro¬ 
ceeds along the following lines. As the ions get heated, 
the dominant wave vector, fc sat , decreases accordingly to 
Eq. (flQ| ) . Since, according to our simulations, (B%) re¬ 
mains essentially constant at late times, there follows 
a steady increase in S p <■ B 2 z )/k 2 sat , and consequently 
of the ion heating. This picture should hold as long as 
the ion anisotropy is sufficient to sustain the instability 
and/or the filament size remains below the transverse size 
of the system. We now propose to derive, from simple 
coalescence arguments [ll],|26j , a closure equation for our 
model describing the temporal evolution of fc sat . 


1. Collective dynamics of the current filaments 

A critical parameter ruling the current filament dy¬ 
namics is the average net current, j x , carried by the 
electron-ion filaments 0. The ion current contribu¬ 
tion, ji Xl can be estimated assuming spatial separation of 
the (initially superposed) counter-streaming ion beams: 
ji x ~ ZieriiVi/ 2, where rii is the total density of the ion 
population (including the two beams). For the reference 
simulation (vi = v e = ±0.2c, rrii = 100ra e ), one thus pre¬ 
dicts ji x ~ 0.1 en c c, in agreement with the current profile 
of Fig. DJd). The electrons tend to neutralize the ion 
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FIG. 9: Screening factor Eq. (l38l) as a function of k sa tc/uj pe 
(black solid line). The estimate j x /jix = 2ck sa t/7ruj pe is su¬ 
perposed as a black dashed line. 


current with an efficiency increasing with the size of the 
filament relative to the electron skin depth. Assuming 
a linear, non-relativistic electron response, the screening 
factor, k = \j x /jix 1 1 of a cylindrical ion current filament 
of diameter A sat /2 = 7r/ fc sat , reads m 


K = 2/i 


/ 7 TU)p e 

V 2c/c sa t 


K i 



(38) 


where I\ and K\ are the modified Bessel functions of 
the first and second kind, respectively. In the limit of 
7ruj pe /2ck sa t 1, the above equation simplifies to 


« - Ji. 


2c/c sat 

7TUJp e 


(39) 


It can be shown that Eq. (139) ) also holds in a planar ge¬ 
ometry. Both expressions are plotted in Fig. [9] Applied 
to the current profiles of Fig. m), exhibiting a typical 
wavelength of A sat cj pe /c ~ 10, the above formula predicts 
a screening electron current of \j ex \ — 0.5 ji x , matching 
the simulation results. From Eq. we deduce the 

approximate net filament current 


jx - ZieriiVi ^^-. (40) 

TTUJpe 

We now derive a differential equation obeyed by the 
average separation length of the filaments, which will be 
be equated to A sat . Furthermore, the diameter of the fil¬ 
aments will be approximated to A sat /2, while their mean 
particle and current density will be assumed constant 
during each merging process. Although similar, our ap¬ 
proach differs from that of Refs, [lj], [26j. The equation 
of motion applied to the distance between two filaments, 
y, reads 


Y ~ 


Jx 


mini 


-(B z ) 


Z{eKVi 

mi 


{Bz), 


(41) 


where (B z ) denotes the averaged magnetic seen by a fil¬ 
ament, which derives from the averaged vector potential 


(A x ). The above equation can then be recast as 

AY^^^A( Ax )(t), (42) 

TTi-i 

1 /2 

Making the approximations e(A x ) ~ S p ' and eA (A x ) ~ 
Sp 1 ^ 2 ASp/2 leads to 


Y 2 ~Z 4 


rS P (t) 
JSM*) 


S P (t*) 2y/5p 


dSp , 


(43) 


where t* denotes the start time of the nonlinear phase 
and we have assumed Y(t*) = 0. The ions are assumed 
to fulfill \£i(t*)\ <C 1. Let us now introduce r c , the typical 
coalescence time between two filaments, such that \Y\ ~ 
Asat/ Fc 0,0. During a merging event, |A sa t| ~ A sa t/r c , 
so that we can estimate |A sat | ~ |y|. There follows 


^sat — %i 


Js p (U ) 


KVi 


*) 


du. 


(44) 


where we have further assumed A sat (£*) = 0. Since this 
equation only involves quantities spatially averaged over 
a large number of filaments, it can be combined to the 
quasilinear equations of Sec. IIIB1 

To make analytical progress, we inject in Eq. (1441) 
simplified forms of Eqs. (I2T]) , ([30]) . ([34]) , (135 ]) and ([38]) , 
in the limits of K 2 ~ rmv § and a* 1 : 


q ~ - / m i v 0 

p “ 2 aiZ? V 2 + ca 


u I 0 

'Gat, — 


1/2 


di — 


2c \ y /2a i Z‘iS P/ 

rriiVo 


Tiy ~ V 0 \jZJZ[s~ p , 

2^2 aiZfSp 


Vi ~ Vo exp 
2c/c sat 


rriiVo 


pi ~ 


TTUOn 


(45) 

(46) 

(47) 

(48) 

(49) 

(50) 


Plugging Eqs. ([45]) - ([50]) into Eq. ([44]) with Vi = vo yields 


A sat — 


UpjV 0 / ,5p(t) du_ ( rUjV 0 

2?r rriiuipe J Sp (t«) V™ \ y/2 a~Zfu 


1/2 


(51) 


The above equation can be readily solved in combination 
with Eq. ([46]) . giving 


^sat — 



where At = t — t* and 


(52) 


TO 


27r(8at) 1/4 / mi 
vo \Zim e 



(53) 
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is the typical coalescence time, that is, the time over 
which fc sat decreases by half. It can also be viewed as the 
lifetime of the slowly-evolving filamentary state estab¬ 
lished at magnetic saturation. As expected, it increases 
with the typical distance between filaments after satura¬ 
tion, A*. Note that the dependence of To upon A* cancels 
out in the long-time limit of A sat (£). Also, the influence 

of the electron screening transpires through the m e 
term. 

The last step consists in substituting Eq. © into 
(|^5|) - (|^9|) to obtain a fully predictive analytical formula¬ 
tion of the plasma parameters as a function of At and of 
the wavevector at the end of the linear phase (fc*): 


4 kic 


2 J2 


CLi ~ 




<4 (l + AtVr 0 2 ) 2 ’ 


pi 

XI 


(1 + A t 2 /r 2 y 


2 aiZf (4fc2 c 2/ w 2. _|_ 2(1 _)_ A 

T ~mv 2 ( 1 + At2 / r o) 2 

w ~ 10 4 k 2 c 2 /^ + 2(1 + AP/tq) 2 ’ 


2 ’ 


Vi ~ Vo exp 


2(1 + At 2 /tq) 2 


4fcJc 2 /Wp i + 2(1 + At 2 /tq) 2 


(54) 

(55) 

(56) 

(57) 


Taking At = 0 gives the plasma parameters at the be¬ 
ginning of the Weibel saturation phase (fc sa t = fc*). 

The cold-limit approximation made in deriving Eqs. 
©-([Sol) (z.e., assuming T^ e (0) = 0) is valid provided 
the initial temperature verifies 


Ti, x , y ( 0 ) < T iy (U) 


rriiVQ 

4^c 2 /^. + 2’ 


(58) 


a condition fulfilled in the previous simulations. 


2. Influence of the initial filament size 


Making use of & = T/fc^ yj2Tio/rrii and S p (t *) ~ (eA*) 2 , 
we obtain 

( 61 ) 

Combining Eqs. m and dal allows us to estimate k* 
as 




UJpi 


k 2 



2aiK 2 

rriiV q ^i 


- 2 . 


(62) 


The term can be estimated by maximizing the growth 
rate computed from the exact dispersion relation (us¬ 
ing the initial plasma parameters). For initially low- 
temperature plasmas, the initial ion anisotropy ratio ver¬ 
ifies af^ rsj rriiVo/T^ 1, so that A* <C c/uj p i. The 
above formulation, based on the Davidson scaling, de¬ 
pends on £^(0) and is thus imprecise. For vi = ±0.2c, 
v e = 0 and rrii = 100m e , linear theory predicts ~ 3.8, 
and hence \*oo pe /c ~ 17 (to be compared with the sim¬ 
ulation result \*uj pe /c ~ 12.6 at uo pe t* = 700) and 
To ~ 4300cj pe . For V{ = ±0.4c and = 1836ra e , one 
obtains £$ ~ 1.8, \*uj pe /c ~ 9.3 with the simulation re¬ 
sult \*LV pe jc ~ 14 at oj pe t* = 1000) and To — lb${)uj pe . 
Yet, the error made in using these estimates should not 
impact the long-term evolution of A sat : 


Asat(At > 3t 0 ) — 2( xi ( ) (iVpiAt ) 2 , 

127r^(2<ai)2 vc/ iv pP 


(63) 


which is independent of A*, yet involves the electron 
mass. Using Eq. m, we can derive the long-time ex¬ 
pression of the ion anisotropy ratio: 


di — 20487r 6 c^ 


rrij 

Zim e 



(64) 


The time required to reach quasi-isotropization (a* = 2) 
can therefore be estimated as 


Our model requires the knowledge of the typical fila¬ 
ment wavelength at the beginning of the nonlinear ion- 
Weibel phase, A*. A crude approximation of A* can be 
made using Davidson’s magnetic trapping model [20j . 
This model assumes that the linear phase of the insta¬ 
bility ceases when the magnetic bounce frequency of the 
driving particles (here the ions) becomes comparable to 
the linear growth rate of the instability. The saturated 
potential vector, A* = A(t *), is therefore expected to 
fulfil 


Tu ~ 


I Zievok^A* 


rrii 


which leads to 


miTj ' 
Zievokl ' 


(59) 


(60) 


Atform 


26__^ / Zjmi \ 1/4 

ujpi vo \ m e ) 


(65) 


where we have assumed a® 1. This time can be 
viewed as a lower limit of the shock formation time since, 
in addition to neglecting the initial electron-driven phase, 
our calculation stops before full isotropization (a* = 0). 
Interestingly, this lower limit scales as mf 4 /my 4 , as a 
result of electron screening. Our formulae will be com¬ 
pared with shock simulations in a forthcoming publica¬ 
tion. 

Figures [TOf a.bi illustrate for two parameter sets the 
theoretical evolution of A sat and cq. In each case, we have 
considered a finite range of values for A*. As expected in 
the di 1 limit, the curves converge to the same limiting 
curve after a few to’s. Finally, Eqs. are plotted 

in Figs. [H16] and [8l where they show overall agreement 
with the corresponding simulation results. 
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(a) 



(b) 



to p. X10 4 

FIG. 10: Temporal evolution of A sa t w ve /c from the numerical 
(red solid line) and analytical approximate (blue solid line) 
resolution of Eq. for two systems: (a) rrii — 100ra e , 

Vi — ±0.2c with the initial conditions 9 < A *uo pe /c < 30 at 
t*uj pe = 600; (b) rm — 1836m e , Vi — ±0.4c with the initial 
conditions 6 < A *u; pe /c < 25 at t*uj pe = 1500. The subpanels 
plot the temporal evolution of the ion anisotropy ratio a*, 
the typical time To, Eq. E3, being indicated by a black 
arrow. The long-time approximations of Eqs. © and (EH 
are superposed as green dashed lines. 

III. CONCLUSIONS 

We have described the slow dynamics of the non¬ 
linear Weibel-filamentation instability using a set of sim¬ 
plified quasilinear relations, valid in the case of highly 
anisotropic ion beams and homogeneous profiles. Fairly 
good agreement between the theoretical expectations and 
2-D/3-D simulations in the non-linear regime has been 
found for various plasma parameters. A closure relation 
modeling the collective filament dynamics has then been 
derived and solved in the high-ion anisotropy limit. Our 


analytical formulae, Eqs. dSl-dSB, are found to cap¬ 
ture with reasonable accuracy the simulation results. We 
have obtained an expression for the “quasi-isotropization 
time” of the ion population, Eq. (EH), which may be 
considered as an upper limit of the shock formation in 
the case of bounded ion beams. An important finding is 
that this time scales as ^“/(m^/me) 1 / 4 due to electron 
screening effects. This result should be taken in con¬ 
sideration when analyzing the results of kinetic simula¬ 
tions run, as is usual, with nonphysical ion [l| or electron 
masses [16|,[l7|. Our non-linear model, albeit based upon 
a number of simplifying assumptions, therefore consti¬ 
tutes a complementary tool to first-principles simulations 
for the understanding of the ion-Weibel-filamentation in¬ 
stability in realistic settings. Its applicability to shock¬ 
relevant configurations will be addressed in a forthcoming 
paper. 
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Appendix A: Electromagnetic dispersion relation for 
bi-Maxwellian distributions 

Let us consider a charge-neutral plasma composed of 
a number of charged particle species (identified by the 
s subscript). The linearized Vlasov-Maxwell equations 
yield the general electromagnetic dispersion relation of 
the plasma between the wave vector, k, and the imagi¬ 
nary frequency, uj (38| : 

e X x ~ k 2 c 2 sin 2 0)(uj 2 £yy — k 2 c 2 cos 2 0) 

— (uj 2 £ xy + k 2 c 2 cos#sin#) 2 = 0 . (Al) 

where 0 denotes the angle between the x-axis and k. The 
dielectric tensor can be expressed in the form 

uj 2 

e aP = Safi + ^ 5 (A2) 

s 

where is the susceptibility tensor and uo ps is the plasma 
frequency of the 5th charged species. In the case of non- 
relativistic bi-Maxwellian distributions [Eq. (|3j)], the ten¬ 
sor elements read 
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Xxx — — 1 + 


cos 1 2 3 4 0 + sin 2 6 ^ 7 ) -sin 20- 

Mi Mx 


M 


/3 2 /iI, + AX sin" 0 — 2 sin 0T- + 2 sin 20-^- 

MX Ml MX_| 2 


[l-£ 2 Z'(0] +2sin^^^-A)-Z(0 
2'(0 


•Za/^A) 


5 0 — sin 0 


M_ 

MX 


tZ'it): 


Xyy ~ 1 

+ 

Xxy = 


M 


cos 2 0—— 2 cos 2 0-^— 2 sin20 

Mx Ml M± 


z\i) 


sin 2 0 + cos 2 O^Ty- + sin 20 — 

Ml 


Mx 

cos 0 sin 0(1 — -^ 2 “) + cos 20— £ [£(£) + £Z'(£)] ——\l ^-/3o cos 0Z(£) 
M_l MxJ Mx V 2 




cos 0 sin 0 


_i_Al 

^ 9 

z'iO 

sin0 + cos 0 — 

Am 


_Mx 

MiJ 

2 

Mx. 

V 


(A3) 


(A4) 


(A5) 


where Z and Z' are the plasma dispersion function and 
its derivative. We have also introduced /3q = vq/c and 


M = cos0sm0 ^ — - — J , 

(A6) 

m cos 2 0 m sin 2 0 

M rri 1 rri » 

-Lx ± y 

(A7) 

m sin 2 0 m cos 2 0 

MX — T T rp 5 

-L x -Ly 

(AS) 

(= \f§(i- v « coal> ) - 

(A9) 

Introducing the normalized wave phase velocity, = 

cj/Mc, Eq. ()A1|) can be recast as 

ak^ T bk 2 T c — 0 , 

(A10) 

with 


a = (/?! — sin 2 0)(/3| — cos 2 0) — cos 2 0 sin 2 0 

, (All) 

b = (sin 2 e - Pi) Z ulsXyy( cos 2 0 - /?!) Z 

s s 

w psXzz 

+ 2 cos 0 sin 0 ^ w p S X 2 /z j 

(A12) 


s 


c = ( Z u lsXyy) (Z U lsXzz) “ ( Z “IsXyz) . 

s s s 

(A13) 


In Eqs. (|A3D - (|A13|h the subscript 8 have been omitted 
on the elements of for the sake of clarity. The wave 
vector is then given by 


,2 _ -KM ± XXX) 
MX) 


(A14) 


with A = vT 2 —4 ac. This formulation, in which the 
squared wave number k 2 {> 0) is a function of /3^ only (for 
a given propagation angle 0), lends itself to the efficient 
numerical scheme introduced by Fried and Gould \39] in a 
non-relativistic electrostatic framework, and generalized 
recently to the electromagnetic regime for various distri¬ 
bution functions [40, l4l| . This scheme consists, first, in 
determining the locus of the zeroes of $sG{P<f>)- This can 
be readily performed by means of a contour plot in a 
finely discretized portion of the complex (3^ plane. Then, 
we retain those zeroes fulfilling 5ftt/(/^) > 0 and identify 
k = y/$lG(/3(f>)- Depending on the /^-domain considered, 
this method allows us to simultaneously solve for a set of 
discrete electromagnetic solutions uo(k,0). 
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